Reptation quantum Monte Carlo for lattice Hamiltonians 
with a directed-update scheme 
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We provide an extension to lattice systems of the reptation quantum Monte Carlo (RQMC) 
algorithm, originally devised for continuous Hamiltonians. For systems affected by the sign problem, 
a method to systematically improve upon the so-called fixed-node approximation is also proposed. 
The generality of the method, which also takes advantage of a canonical worm algorithm scheme 
to measure off-diagonal observables, makes it applicable to a vast variety of quantum systems and 
eases the study of their ground-state and excited-states properties. As a case study, we investigate 
the quantum dynamics of the one-dimensional Heisenberg model and we provide accurate estimates 
of the ground-state energy of the two-dimensional fermionic Ifubbard model. 



I. INTRODUCTION 

The path-integral formulation of quantum mechanics 
is the foundation of many numerical methods that allow 
one to study with great accuracy the rich physics of inter- 
acting quantum systems. At finite temperature, a path- 
integral Monte Carlo (PIMC) technique for continuous 
systems has been developed and applied by Ceperley and 
PoUockji'Si Recently, this approach has been renovated in 
a new class of methods known as worm algorithms'^ 
The zero-temperature counterparts of the PIMC algo- 
rithm are the reptation quantum Monte Carlo (RQMC)^ 
and the path-integral ground state (PIGS) methods,^ 
which have been demonstrated useful to simulate cou- 
pled electron-ion systems^S as well as to infer spectral 
properties from imaginary time dynamics.'* A number of 
important physical problems -particularly in the fields 
of strongly correlated fermions and cold atoms- can be 
fruitfully modeled by lattice Hamiltonians. A first appli- 
cation of path-integral techniques to (boson) lattice mod- 
els was proposed by Krauth et al. in 1991.— Few other at- 
tempts to apply PIMC to lattice models have been made 
ever since, with a recent application of the RQMC idea to 
the quantum dimer model Hamiltonian.^*' In this paper, 
we propose a new method that generalizes and improves 
the approach of Ref. [l^ in several ways. Our method 
is based on continuous-time random walks and is there- 
fore unaffected by time-step errors. Inspired by the work 
of Syljuasen and Sandvikii and Rousseau^^ we adopt a 
generalization of the bounce algorithm of Pierleoni and 
Ceperley,~ called directed updates, which helps reducing 
the correlation time in path sampling. We also intro- 
duce a worm-algorithm based method to calculate pure 
expectation values of arbitrary off-diagonal observables, 
which are generally out of the scope of existing lattice 
ground-state methods. 

The resulting algorithm naturally applies to fermions, 
using the fixed-node approximation. A technique to 
improve systematically upon this approximation is pro- 
posed, based on the calculation of a few moments of the 



Hamiltonian. Our methodology is demonstrated by a few 
case studies on the one-dimensional Heisenberg and the 
two-dimensional fermion Hubbard models. 

This paper is organized as follow: in SeclUwe present 
the general formalism of ground-state PIMC for lattice 
models; in Sec. Illll our implementation of the RQMC al- 
gorithm on a lattice is presented. In particular, we give a 
detailed account of the above mentioned directed update 
technique (Sec. IIII Ap and of the continuous-time propa- 
gator fSec. lIII B|) : in Sec. lIII Cl we introduce an extension 
of the algorithm to cope with off-diagonal observables, 
while in Sec. lIIIDl a further extension to systems affected 
by sign problems is presented, including a strategy to 
improve systematically upon the fixed-node approxima- 
tion. Sec. IIVI contains a few case applications, including 
the simulation of the spectral properties and spin corre- 
lations of the one-dimensional Heisenberg model and the 
calculation of the ground-state energies of the fermionic 
Hubbard model with a significantly better accuracy than 
that achieved by the fixed- node approximation. In Sec.lVl 
we finally draw our conclusions. 



II. GENERAL FORMALISM 

Let us consider a generic lattice Hamiltonian H and 
a complete and orthogonal basis set, whose states are 
denoted by \x). Given the generic wave function \^), 
its amplitude on the configuration \x) will be denoted 
by ^'(x), namely ^'(x) — {x\^). The exact ground-state 
wave function j^I'o) can be obtained by the imaginary 
time evolution of a given variational state |5'y): 



l*o) 



cx lim 



(1) 



where l^*^) = e^'^^l^y), provided that the variational 
state is non-orthogonal to I'I'o), i.e., {'^v\^o) 7^ 0. Then, 
the ground-state expectation value of a quantum opera- 
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tor O can be obtained by 



(O) = lini 



(2) 



A practical computational scheme can be conveniently 
introduced by considering a path-integral representation 
of the imaginary time evolution. To such a purpose, we 
split the total imaginary time /3 into M slices of "du- 
ration" T = /3/Af, in such a way that the value of the 
evolved wave function on a generic many-body state of 
the system reads 



M 



(3) 



XI...XM i — 1 



where we have introduced the imaginary time propaga- 
tors 



(71 



{Xi-i\e 



(4) 



Within this approach, it is easy to write expectation val- 
ues of operators O that are diagonal in the chosen basis 
jx), i.e., {x\0\y) — 0(x)5x,y In fact, in this case we have 
that: 



Exn^(x)0(xM) 



(5) 



where the summation is extended to all possible imag- 
inary time paths X = {xq.xi, . . . ,X2Af}, and 11^ (X) is 
given by: 



n^(x) = ^v{xo) 



2M 



'i>v{x2M)- (6) 



The ground-state energy can be conveniently obtained by 
means of the mixed average: 



Eq — lim 



(7) 



where El{x) = |\I'v')/(a;|^v) is the so-called local 
energy. 

Besides the static (i.e., equal-time) correlation func- 
tions, this formalism allows one to calculate also dy- 
namical correlations in imaginary time Cab{T) = 
{A{T)B{Q)) that can be computed as 



Cab{T) 



y Exn^(X)A(x„)j3(a;„) 



where Xn and Xm are two coordinates of the path such 
that (m — njT = T . 



III. REPTATION QUANTUM MONTE CARLO 



whenever 11^ (X) > for all configurations X. Indeed, in 
this case, 11^ (X) can be interpreted as a probability dis- 
tribution that may be readily sampled by using Monte 
Carlo algorithms. This fact allows ground-state expec- 
tation values to be calculated exactly, within statistical 
errors. 

The basic idea of the RQMC algorithm is to sample 
the distribution probability 11^ (X) by using a Markov 
process with simple moves. Given the configuration X = 
{xo, a;i , . . . , X2m}, a new configuration is proposed in two 
possible ways: either Xl = {xt,xq, . . . ,X2M-i} (which 
we call "left move") or X^j = {xi, . . . , X2m, xt}, (which 
we call "right move"). In both cases, xt is a new config- 
uration proposed according to a suitable transition prob- 
ability {x — Xt), where x stays for xq {x2m) when 
the left (right) move is considered. Such "sliding moves" 
are depicted in Fig. [1] Ideally, the transition probability 
should guarantee the minimum possible statistical error 
on the desired observables and, to such a purpose, it has 
been proved useful to consider the propagator with im- 
portance samphng, i.e., G^y = G^.y'^ v {y) / v (x) and 
take the following transition probability 



where 



R^{x 



w{x) 



w[x) 



(9) 



(10) 



represents the normalization factor. The explicit form of 
R'^{x — ^ Xt) will be discussed in more detail in Sec. lIIIBl 
The proposed configuration X^; (where d = L or R) is 
accepted or rejected according to the usual Metropolis 
algorithm, where the acceptance rate is given by: 



A ~ min < 1 



Il^(Kd)R-ixT^x) 
' nP{X)R-^{x ^ xt) 



(11) 



In this way, a sequence of configurations X'^ is generated, 
k being the (discrete) time index of the Markov chain. 

In order to reduce the auto-correlation time of the ob- 
servables it is convenient to make several consecutive slid- 
ing moves along the same imaginary time direction.^ To 
such a purpose, a recent development called "bounce" al- 
gorithm has been proposed<i Although the bounce algo- 
rithm sampling procedure does not fulfill the microscopic 
detailed balance, the equilibrium probability n'^(X) is 
correctly sampledi^ The RQMC algorithm with bounce 
moves can be then summarized in the following steps: 

1. For the current direction of the move and for the 
present configuration X'^ , propose xt according to 
the transition probability R'^{x — > xt), where x = 



Xq ii d = L and x 



if d = i?. 



A probabilistic interpretation of the previous expecta- 
tion values ([5]), (O, and ([S]) can be immediately recovered 



2. Given the form of the acceptance ratio A of 
Eq. (ITT|) . accept the proposed configuration accord- 
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Figure 1: Pictorial representation of the "sliding moves" along 
the right imaginary time direction. In the new configuration 
(bottom), a new head for the reptile is generated from the old 
configuration (top) and the tail is discarded. 



ing to the probability 



A, 



min < 1 



\i d = L, or with probability 



A 



R = mm 



(12) 



(13) 



3. If the move is accepted, update the path config- 
urations according to X*^"*"^ = and continue 
along the same direction, otherwise X^-'^^ = X*"' 
and change direction. 

4. Go back to 1. 



A. Directed updates 

At this point we introduce a novel alternative sam- 
pling approach, which generalizes the bounce idea while 
strictly fulfilling the detailed balance condition. Such a 
scheme, which is largely inspired by the loop algorithm 
methods devised for the stochastic series expansio n^^i^^ 
and for the worm algorithmf^^^ allows one to choose the 
time direction in a purely Markovian way, i.e., indepen- 
dently of the previous history. 

In our algorithm, a Markov step consists of many sim- 
ple consecutive "sliding moves", whose number is not 
fixed a-priori but is determined by a certain probabil- 
ity (see below). The actual Monte Carlo step takes place 
at the end of few consecutive updates along the currently 
chosen direction. In the examples below, we denote the 
number of these sliding moves between two Monte Carlo 
steps by s. 



At the beginning of each Markov step we choose a di- 
rection d according to the probability P(X'^,d), whose 
form will be specified later. Assuming that the right di- 
rection has been chosen, we propose a new configuration 
XT, according to the transition probability R'^{x2j^j — >■ 
xt) and the configuration labels are shifted according 



^2M 



XT- 



At 



to X'^+i = {x\,. . . ^x^mi^^t], with 
this point, we continue updates along this direction 
with probability i4r(X'^+^, — or stop with probability 
[1 — if (X'^^^, —?►)]. If it has been decided to continue the 
updates, then a new configuration is generated according 
to i?'^(a;2M ~^ ^t) and the labels of the configuration are 
again shifted according to X'^+^ = {xj^"'"^, . . . , x\'l^ ,xt}- 
The Markov step finishes after s consecutive updates 
along the right direction only when X(X'^+^, — >■) < 
where is a random number uniformly distributed in 
[0, 1). At this point a Metropolis test should be done, in 
order to accept or reject the sequence of intermediate s 
sliding moves: 



A = min < 1, 



where (see Appendix E| 



q(X'^+^) 



(14) 



9(X) 



^(x,^) 

l-if(X,^) 

^(x,^) 

l-if(X,^) 



w{x2M-i) 
w{xi). 



(15) 



However, in order to avoid time-consuming restorations 
of the original configuration, it is preferable to accept 
all the moves, while keeping track of the residual weight 
(7(X). This is possible since A only depends upon ini- 
tial and final configurations, so that, given that all the 
intermediate moves are accepted, the sampled distribu- 
tion probability is 11^ (X) x g(X). The contribution of 
the current configuration to statistical averages must be 
then weighted by the factor l/q{yi). 

To proceed to the next Markov step, a new direction d 
is chosen according to P(X'^+'*, ^) and P(X'''+'', — ^) and 
the updates are carried along the extracted new direction. 

Let us now show the actual expressions for the afore- 
mentioned probabilities. In Appendix |21 it is demon- 
strated that the detailed balance is satisfied if one chooses 
the probabilities for the directions as 



p(x,H = 

P(X,^) = 



1 



1 + a(X) ' 
«(X) 



where 



a(X) 



1 + a(X) ' 

w(x2M~i)l~K{y.,^) 
w{xi) 1-X(X,^)' 



(16) 
(17) 

(18) 



which is positive and, therefore, guarantees that the 
above defined quantities are well defined probabilities. 
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i.e., < P(X,^) < 1 and < P(X, ^) < 1, with the 
additional property that P(X, ^) + P(X, = 1. 

Regarding the probabihties to continue the updates 
along the current direction, we have a substantial free- 
dom of choice, provided that the condition = 

w(ji>2m\) ' satisfied (see Appendix E]) . In this paper we 
have adopted 



i^(X,^) 

i4r(X, — >) — amin-^l. 



a min {1, fe(X)} , 
1 



6(X) 



where we have defined: 



6(X) 



w{x2M~l) ' 



(19) 
(20) 

(21) 



and < Of < 1 is an arbitrary parameter of the algo- 
rithm, which controls the average number of consecutive 
updates along the same direction. 

Summarizing, the RQMC algorithm with directed up- 
dates consists of a sequence of Markov steps determined 
by the following rules: 

1. Choose a time direction d according to the proba- 
bilities of Eqs. (dni) and P71) . 

2. Propose a new configuration xt according to the 
transition probability K^i^x — >■ Xt), where X — 

a d = L and x = if = -R- 

3. Shift the configuration indexes according to 
X'=+i = {xt, xg, . . . , x^M-il a d = L or X'^+i = 



{x1 



Xt} if d = i?. 



4. According to the probability K{'X'',^) or 
iiT (X*^, -i— ), decide whether keep moving in the same 
direction or change direction. In the former case, 
go to 2, otherwise go to 5. 

5. The Markov step ends here and the current configu- 
ration carries the weight l/q{'X.''^^), where s is the 
number of intermediate moves along the direction 
chosen. 

The relationship between the directed update scheme 
and the bounce algorithm is further elucidated in Ap- 
pendix [B] where general considerations about the effi- 
ciency of the algorithms are presented. 



B. Continuous-time propagator 

One of the most striking differences between the origi- 
nal formulation of the RQMC on the continuum and the 
present formulation on the lattice is the lack of the dis- 
cretization error appearing in the Trotter decomposition 
of the propagator. Indeed it is easier to carry the propa- 
gation in continuous imaginary time on a lattice,J^ than 
on the continuum.^® To such a purpose, let us consider 



the limit of an infinitesimal imaginary time e, for which 
the transition probability of Eq. ^ can be written as 



R'ix^y) 



5xy [l + eEL{x)]-e 



eEL{x) 



H 



xy 



^'v{x) 



(22) 



oie'), (23) 



where El{x) is the previously defined local energy and 
Hx,y = {x\H\y) denotes the matrix element of the Hamil- 
tonian. Whenever "^v{y)Hxy/'^v{x) is non positive for 
all X and y, this equation takes the form of a continuous- 
time Markov process, whose analytical properties are well 
known. In particular, the probability distribution for the 
"waiting time" in a given state x, i.e., the average time 
that the system spends in the state x before making an 
off-diagonal transition to another state y ^ x^ is exactly 
known, namely P{tw;x) = exp{-r^ [Hxx - El{x)]}. As 
a consequence, the finite-time propagator {x — > y) can 
be directly sampled, giving rise to a succession of a cer- 
tain number n of consecutive transitions a: — > zi — >■ Z2 — >■ 
•••—>■?;, with corresponding waiting times r^(zi) (such 
that 'Ylii'''w{zi) — t). The normalization of the whole 
process is 



w{x) — cxp 



,{zi)EL{zi) 



(24) 



where the waiting times are extracted according to 
the exponential probability P{Tw',Zi). The transi- 
tions between the intermediate configurations are done 
according to the off-diagonal elements of Eq. (1^51) . 
i.e., Zi+i is chosen with probability proportional to 



C. Off-diagonal observables 

The formalism so-far developed allows one to success- 
fully compute pure ground-state expectation values of op- 
erators that are diagonal in the local basis x, with the 
expectation values of off-diagonal operators restricted to 
the so-called mixed averages^'^^ Nonetheless, it is of- 
ten of great interest to remove such a limitation (whose 
result is biased by the quality of the variational wave 
function) and a dedicated sampling strategy has to be 
devised in order to cope with such a need. In the follow- 
ing, we show that a relatively easy modification of the 
sampling scheme can accomplish this task, providing us 
with a general tool to compute ground-state averages of 
operators that are non local in the chosen basis x. 

Let us consider an arbitrary off-diagonal observable O, 
which can be in turn considered as the summation of 
many observables we are interested in, i.e., O = 0'^'^\ 
For example, we can imagine these operators to be the 
components of the one-body density matrix at a given 
distance, — X](r r') ^r^"^' ^'^^^ the summation ex- 

tended to all lattice coordinates at a fixed distance d. 



In the spirit of Refs |12||1 
operator defined by 



we introduce a worm- 



(25) 



where A is a positive constant, and consider the extended 
configuration space spanned by the paths 



4=1 

2M+1 

X n G:^_^x^x^v{x2m+i). (26) 

i=R+l 

The extended paths are broken in two (space)- 
discontinuous pieces by the worm operator, which is 
placed at an imaginary-time < t^h < /?. Therefore, 
paths contain 2(M + 1) coordinates, including xl and 
xji that refer to the same imaginary time r^fl. 

The configuration space spanned by Eq. is clearly 
larger than the one spanned by Eq. ([6]), which is recov- 
ered whenever xl = xr, i.e., when the worm operator is 
diagonal. 

The pure ground-state expectation value of the oper- 
ator O is conveniently written in terms of the extended 
paths as 



0. ^ i Exn^(x)xe(x^^x^) 

\ 13- 



Exn(^(x)xe(.Ti 



xr) 



(27) 



where 6(C) 7^ whenever condition C is satisfied. The 
modulus of Eq. (1261) can be in turn interpreted as a prob- 
ability distribution and stochastically sampled by means 
of the elementary sliding moves considered before. In- 
deed, whenever the worm operator is far from the ends of 
the imaginary-time paths, the sampling scheme remains 
unchanged. In this case, a move along direction d will 
generate a new head (or tail) for the reptile according to 
W{x — >■ xt) while shifting the worm position of ±t. On 
the other hand, whenever the worm operator reaches the 
ends of the reptile, a new worm configuration is proposed 
on the other side; in analogy with the previous analysis, 
new configurations are generated according to a transi- 
tion probability 



R'^ix ^ y) 



1 



■w{x) 



V'v(y) 



ipv{x) 



(28) 



where w(x) is the normalization factor. Due to the 
particular form of the matrix elements (j25p . the tran- 
sition probability will lead either to diagonal config- 
urations {x = y) or to off-diagonal configurations 
{x 7^ y), thus generating continuous and discontinu- 
ous paths. The relative probability for diagonal and 
off-diagonal configurations depends on the value of A 
that can be tuned in order to reach a balanced sam- 
pling frequency for the difi'erent sectors of the extended 
paths. In order to exemplify the worm updates, let 



X 




Xl 



X2 



Xi 



X4 



x^ 




Xl 



X2 



XZ 



Xii 



XT 



Figure 2: Pictorial representation of the "sliding moves" along 
the right imaginary-time direction when the worm operator 
sits at the tail of the reptile. In the new configuration (bot- 
tom), a new head for the reptile is generated from the old 
configuration (top), the old tail configuration is discarded and 
the worm discontinuity is moved to the "neck" of the reptile. 



us consider the case in which d — R and a configu- 



ration '^v [xo) yVx 



I-2M-I-1 



"^v {x2M-\-i), af- 



ni=2 ^ Xi-ix 

ter a sliding update in the right direction, we will have 

2Af+i2:T^v [xt), where xt 
is proposed according to the transition probability 
-R^(a;2A/+i Xt) (see Fig. [2]). In analogy with the pre- 
vious case, the acceptance factor for the bounce moves 

reads ^«^min{l,fef)}. 

Summarizing, the RQMC with worm-updates consists 
of the following steps: 

1. For the current direction of the move d and for 
the present configuration X*^ consider the worm- 
operator position tl/j. 

2. If the worm is not at the ends of the reptile (i.e., 
tlr 7^ when d — L and ttr 7^ /? when d = R) go 
to step (a), otherwise go to step (b). 

(a) Propose a new configuration xt according to 
the transition probability R'^{x — >■ cct), where 
a; = ccq if d = L and x — xl^jy^^^ ii d = R. The 
new configuration is accepted with probability 

w{x^) 



At = min < 1 , 



w{x2m) 



if d = L, or with probability 



min < 1 



^(^2M+l) 

w{x'[) 



(29) 



(30) 



if d = i?. In the proposed state X^, all the 
configuration labels are shifted in the d direc- 
tion, determining in turn a shift of the worm 
operator of a time interval ±t, depending on 
d. 
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(b) Propose a new configuration xt according to 
the worm transition probability FO^ [x —^xt), 



where x 



ii d = L and x 



^2M+i d = R. Accept the new configuration 
with probabihty 



Al = min <! 1 , 



wix^M) 



if d = L, or with probabihty 



A]^ = min < f , 



''^(^2M+l) 

w{xi) 



(31) 



(32) 



if d = i?. In the proposed state X^, ah the 
configuration labels are shifted in the d direc- 
tion, and the worm operator is moved from 
the head (tail) to the tail (head) of the rep- 
tile, depending on d. 

3. If the move is accepted, update the path config- 
urations according to X'^^^ = X^ and continue 
along the same direction, otherwise X^-'^^ — X*"' 
and change direction. 

4. Go back to 1. 

This scheme samples the probability density associated 
to the modulus of Eq. (pS)) . and the expectation values 
of the individual components Od can be recast as statis- 
tical averages over such a probability distribution, while 
keeping track of the overall sign of the extended paths. 
In particular the best estimate of the ground-state expec- 
tation values is obtained when the worm is in the central 
part of the path, at tlr — 13/2, leading to 



Exn^(x) X e (oilU ^ 0, = f 
Sxnw(^) X e (^XL = XR, = 



^ (e (o^Sn ^0]x sign 
A 



n^(x) 



center 
OD 



/\reentei' 



(33) 



where (. . . )q^'^^ denotes statistical averages over the off- 
diagonal distribution 



and ]V|j°"'°'' is the number of configurations sampled with 
a diagonal worm operator in the center of the paths. 



D. Tackling the sign problem 

When the probability distribution of Eq. © is not pos- 
itive defined, as is generally the case with fermions, the 
probabilistic interpretation of the imaginary time paths 
breaks down. This circumstance, which is known as "sign 
problem", originates whenever 'i/v{y)Hxy/'i'v{x) > for 
some element x ^ y. In this case, it is not possible 
to have polynomial algorithms that are able to obtain 



an exact solution of the problem, which would imply to 
sample correctly the resulting signs. Therefore, approxi- 
mated schemes are welcome and often adopted, the most 
widespread one being the so-called fixed-node (FN) ap- 
proximation. For lattice systems, this approach relies on 
the definition of an effective Hamiltonian, which depends 
parametrically on the nodal structure of a given vari- 
ational wave function 5'y(a;) — (a;|5'y)iA- The matrix 
elements of the FN Hamiltonian are defined as 



H 



in 



Hx 




xy 



+ Vsi{x) if X = y 

ii^v{y)Hxy^v{x)<0 (34) 
if *v(y)iJ:ry*y(x) > 



where the sign-fiip potential is Vsi{x) = 
J2y:sf'^viy)Hxy/^v{x), the sum being extended 
to all the sign-flip states defined by the condition 
'ifv{y)Hxy^v{x) > 0. With such a choice, the transition 
matrix i?' (x — )■ y) of Eq. (j23|) is always positive definite 
and the summation of Eq. Q is now restricted -which 
results in the FN approximation- to a region of the 
Hilbert space in which imaginary time paths are positive 
definite. Therefore, within the FN approximation, the 
ground-state wave function of W'^ can be stochas- 

tically sampled without any sign problem. Moreover, 
it is easy to show that the FN approximation becomes 
exact whenever the signs of the exact ground state are 
known. Most importantly, it has been provenii that the 
FN ground-state energy E'^" = (i?^") gives a rigorous 
upper-bound to the exact ground-state one and improves 
the pure variational results. 

At this point, we introduce a straightforward, although 
computationally expensive, way to improve further the 
FN energy. Our strategy amounts to compute the expec- 
tation values of arbitrary powers of the original Hamilto- 
nian H on the FN ground state 1 , namely 



*fn|i?'=l* 



(^'f„|*fn) 



(35) 



The FN ground state can be expanded in the basis set of 
the eigenstates of iJ as l^-fn) = 7o|*o)+7i|^'i)+72|^'2) + 
. . . and Lfc = 7o^^o +7?^f +7|£^l + ■ • ■ , with Z^ ll = ^- 
Since very often the FN wave function has a con- 
siderable overlap with only few low-energy states, the 
knowledge of the first few moments of the Hamiltonian 
are enough to approximately reconstruct both the coef- 
ficients 7i and the energies Ei. To such a purpose, let 
us consider a typical situation in which only the first 2n 
moments of the Hamiltonian have been numerically cal- 
culated and are therefore known. We can then truncate 
the expansion for to the order n — 1 having a closed 
system of 2n equations 



i=0 



■7^ 

n.n i.n J 



(36) 



for fc = 0, . . . 2n — 1 that can be solved for the unknowns 
7i^„ and In the limit of large n, the approximated 
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i?o,n converges to the exact ground-state energy. More- 
over, we verified that i?o,n > Eq, as a result of a con- 
nection between the solutions of the Eq. (l36l) and the 
Lanczos procedure written in terms of the moments of 
the Hamiltonian.^^ 

The Hamiltonian moments are off-diagonal operators 
and can, in principle, measured according to the sam- 
pling procedure detailed in Sec. IIII CI In the present im- 
plementation we are able to achieve sufficient statistical 
accuracy only for the first moment of the Hamiltonian, 
i.e., Li — (-ff), while higher moments are too noisy. Yet, 
to our knowledge our algorithm is the only one that al- 
lows the calculation of the expectation value of the origi- 
nal Hamiltonian H. This is knownii to be a better upper 
bound than the expectation value of the FN Hamiltonian 
accessible with other zero-temperature algorithms. 

Although we are not currently in position to measure 
directly the Hamiltonian moments Lk we have a con- 
trolled access to the mixed averages 



(*fn|*V') 



(37) 



which present optimal statistical uncertainty. More- 
over, an improved estimate of the ground-state energy 
based on the knowledge of the first few moments i™"^ 
can be obtained solving a system of equation similar to 
Eq. (1551) that leads to the approximate ground-state en- 
ergies E™^^. Unfortunately, the proof that E™^ > Eq, 
for n > 1, is far from being trivial, requiring a general- 
ization of the already non-trivial upper bound for n = 1 
described in Ref. Nonetheless, we have numerically 
verified that, in all the cases treated in this paper (where 
Eq is a-priori known) , the condition E'™i^ > Eq is always 
verified. We are then led to conjecture that this may 
always be the case. 



IV. RESULTS 



<1 




Figure 3: Lowest-energy excitations as a function of the wave- 
vector q for an L = 20 Heisenberg chain. The energies are 
extracted from the dynamical structure factor S{q,ij) and are 
compared to exact results by the Lanczos method. 




Figure 4: The same as in Fig. [3] for L ■ 
given by Bethe ansatz. 



80. Exact results are 



A. Low-energy excitations and spin correlations of 
the Heisenberg model 

Hereafter, we present a simple application of the previ- 
ous ideas to sign-problem free spin Hamiltonians. Let us 
consider the one-dimensional quantum Heisenberg model 

7^ = J^S,-S,+i, (38) 

i 

where = (^Sf, S'f , S'f ^ is the spin 1/2 operator on the 

site i and J > is the nearest-neighbor super-exchange 
coupling. 

The total number of sites is denoted by L and periodic- 
boundary conditions are assumed. This model can be 
solved exactly by using the so-called Bethe ansatz tech- 
nique.^^ Information on the excitation spectrum can be 



obtained from the dynamical structure factor 

5(<z,^) = I dt{Sl{t)Sl^me^-\ (39) 

where S^{t) — l/vT^^ S'|(t)e*'^^ is the Fourier trans- 
form of time-evolved spin projection on the z-axis. By 
introducing a complete set of eigenstates of the Hamilto- 
nian \^n) with eigenvalues E^ we have that 

S{q,Uj) = J2 |(*0|^||*n)P^(^ - C^n), (40) 

where ujn — {E„ — Eq). In the thermodynamic limit, 
the spin-1 states form a branch, which is very similar 
to spin waves in standard ordered systems, although no 
long-range order is found in one dimension. 

Imaginary time correlation functions of arbitrary (di- 
agonal) operators can be efficiently evaluated via Eq. ([S]) . 
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2. 
?3' 




15 20 25 30 35 40 
d 



Figure 5: Ground-state expectation value of the spin-spin cor- 
relation function C(d) for the Heisenberg model on a 80-site 
chain. 



In Fig. [51 we show the expectation value of C{d) for a 
80-site one-dimensional lattice. In this case we are able 
to achieve very good statistics for the off-diagonal ob- 
servable, with a relatively negligible computational effort, 
when compared to the evaluation of the ground-state ex- 
pectation value of other diagonal observables. 



B. Ground-state properties of the fermionic 
Hubbard model 



As an example of the application of the RQMC to sign- 
problem affected Hamiltonians, we present some results 
for the fermionic Hubbard model on a square lattice, de- 
fined by: 

= -i ^ 4<tCj\<t + /i-c. + C/^n^,tf^^,;, (43) 



This fact allows us to have a direct access to S{q, T) = 
{Sg{T)Stq{0)). This imaginary time correlation func- 
tion can be then analytically continued, by using the 
Maximum-Entropy method,^- in order to have a reason- 
able numerical estimate for the dynamical structure fac- 
tor of Eq. (gOl). 

Before presenting the results, let us mention that we 
consider the following Jastrow state as a variational wave 
function)^i2^ 



exp 



\FM) 



(41) 



where \FM) is the ferromagnetic state along the x direc- 
tion, for which {x\FM) does not depend upon the spin 
configuration and the variational parameters Vij are op- 
timized by using the method of Ref. 23. 

In Fig. [21 we show the results for a small i = 20 sys- 
tem, where exact diagonalizations are possible by using 
the Lanczos method. We report the energy excitations 
AE{q) — Eg ~ Eq for the lowest state with 5 = 1 and 
fixed momentum q. In this case a perfect agreement be- 
tween our RQMC results and the exact ones is found. 
Moreover, also on larger systems a very good accuracy is 
possible (see Fig. [4|) , demonstrating the performances of 
our numerical algorithm. 

In order to exemplify the potentialities of the scheme 
outlined in llll C[ we conclude this part of the results de- 
voted to the Heisenberg model showing the ground-state 
expectation value of the spin-spin correlation at distance 
d 



C{d) 



S ,; • S ,; 



(42) 



The desired observable is used as a worm operator and 
the value of the correlation function at the various dis- 
tances is computed by means of the estimator of Eq. 



where (...) indicate nearest-neighbor sites, cj ^ (ci^o-) cre- 
ates (destroys) an electron on the site i with spin tr, and 



at 



As a variational state we consider 



\'^v) = cxp 



\FS) 



(44) 



where \FS) is the non-interacting Fermi sea and the Jas- 
trow factor involves density-density correlations. The 
variational parameters Vij entering in the Jastrow factor 
may be optimized again by minimizing the variational 
energy with the method of Ref. [2^. In order to avoid 
open shells in \FS), we consider 45-degrees tilted lattices 
with L = 2 X sites, such that both the half-filled case 
and selected holes-doped cases are closed shells. 



0.05 
0.04 



^ 0.02 
<1 

0.01 I- 




Figure 6: Ground-state energy for the fermionic Hubbard 
model at half filling on a 18-sites tilted-square lattice. The 
energy difference AE — -Eexact — -E is computed with distinct 
approximations described in the text. 
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U/t 




{^^> 




4 


-42.850(1) 


-43.16(1) 


-43.282(1) 


5 


-36.364(1) 


-36.51(1) 


-37.052(1) 


6 


-31.885(1) 


-32.17(1) 


-32.640(1) 


7 


-28.318(1) 


-28.66(1) 


-29.022(1) 


8 


-25.382(1) 


-25.62(1) 


-26.056(1) 



Table I: Ground-state energy as a function of the Hubbard U 
repulsion on the 50-site lattice at half filling. 



N 



rpmix 



Eai 



50 -42.850(1) -43.16(1) -43.282(1) -43.983(1) 

42 -53.402(1) -53.57(1) -53.769(1) -54.001(1) 

26 -55.4325(1) -55.63(1) -55.6112(1) -55.782(1) 

18 -50.4127(1) -50.50(1) -50.4383(1) -50.474(1) 

Table II: Ground-state energy as a function of the number of 
electrons A*' for Hubbard repulsion C//t = 4 on a 50-site lattice. 
The numerically exact results obtained by the Auxiliary-Field 
Monte Carlo method Eaf are also shown for comparison. 



Let us start by showing the resuhs for 18 electrons 
on 18 sites, where Lanczos diagonalizations are possi- 
blei^ In Fig. ^ we report our results for the ground- 
state energy. The FN approach gives rather accurate 
results for small values of U/t, i.e., U/t < 4, where 
(£^cxact - £^^")/£^exact < 0.01. By Increasing the on-site 
interaction, the FN approach becomes worse and worse. 
This fact is due to the choice of the variational wave func- 
tion that does not contain antiferromagnetic order. Re- 
markably, a considerable improvement may be obtained 
by considering the pure expectation value of the Hamilto- 
nian, which is systematically lower than the FN energy, as 
demonstrated in Ref . [l3 and now accessible within our al- 
gorithm. Further improvements to the FN energy can be 
obtained upon considering few (up to three) higher mo- 
ments of the Hamiltonian measured as mixed-averages, 
see Fig. [S] The scheme based upon the Hamiltonian mo- 
ments (described in Sec. IIII D|) allows us to reach a great 
accuracy for the ground state energy, with a residual er- 
ror almost independent of U/t. Indeed, in this way we 
have (£;cxact - £^)/£^cxact < 0.002 up to U/t = 8. 

This approach remains very effective also for larger sys- 
tems, even though the variational wave function loses ac- 
curacy by increasing the cluster size (because the ground 
state has antiferromagnetic order in the thermodynamic 
limit, while the variational state is paramagnetic). In 
Table HI we report the ground-state energy for 50 sites 
for the half- filled case, while in Table HIl we report the 
ground-state energies for selected cases at finite hole- 
doping, where numerically exact results (for moderate 
values of U and moderate lattice sizes) can be obtained 
by the Auxiliary-Field Monte Carlo method/^ 



V. CONCLUSIONS 

In this paper we have provided an efficient and gen- 
eral formulation of the reptation quantum Monte Carlo 
technique on lattice models. In particular, we showed 
an alternative sampling approach which generalizes the 
bounce algorithm, previously introduced to reduce auto- 
correlation time of the observables. Our scheme allows 
one to choose the time direction in a purely Markovian 
way. In addition, the average number of consecutive 
moves along the time directions may be optimized by 
a fine tuning of a certain parameter that has been ex- 
pressly introduced in the transition probabilities. We re- 
ported benchmarks for two different models with pure 
bosonic and fermionic degrees of freedom, by showing to 
what extent it is possible to have accurate results both on 
the ground state and low-energy excitations. The intro- 
duction of a general method to compute ground-state ex- 
pectation values of arbitrary off-diagonal observables also 
constitutes an important achievement, which will ease the 
study of relevant properties such as Bose-Einstein con- 
densation and superconductivity phenomena in strongly 
interacting models. In addition, the possibility to di- 
rectly measure the pure ground-state expectation values 
may open the way to a better optimization of the corre- 
lated wave function associated to the ground-state of an 
effective Hamiltonian which is not the FN one. 
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Appendix A: Derivation of the probabilities for the 
directed-update scheme 

In this Appendix we give a detailed derivation of the 
probabilities for the directed updates. The detailed bal- 
ance condition guarantees that the given probability dis- 
tribution H'^(X) is sampled if transitions from an initial 
state X*"' to a final state X*"'"*"" differing for s intermediate 
updates are accepted according to: 



A' 



min < 1 



H'3(x'=+") r"(x''-+'^ ^ X* 



(Al) 



being the overall transition probability between the 
two states. Let us first consider the case when s = 1 
and fix the right direction d = R {a, similar derivation 
can be obtained for d = L). In this case, the transition 
probability from the initial state to the final state reads 



r^(x'^ 



p(x^^)xi^-(x2\, 

[l-X(X^+i,^)], 



^2M ) 



(A2) 



10 



namely, it is the product of the probabihty of having cho- 
sen the right direction, times the transition probabihty 
for the new tail of the reptile, times the probability of 
stopping the updates after one intermediate step. The 
inverse transition probability instead reads 



k+l 



X 



= P(X'^-+1,^) X ^ x'^) X 

X [l-i^(x^4-)], (A3) 



which can be obtained reversing the time directions and 
considering transitions from the head of the reptile in- 
stead that from the tail. Therefore, the acceptance factor 
reads as 

. 1- jr(x^H 

mm < 1, — ^ — X 



P(Xfc,^) 
fc+i ^ p(x''^+i^ 



^T^ — 



1 - ^'(xfe+i,^) 



(A4) 



For two intermediate transitions instead the transition 
probabilities are 



r^(x'^ ^ x" 



fc+2\ _ 



= P(X^^) X 



2M 



^2M ) ^ 



X if(X^-+\^)xi?-(x^+/^x2^+/)x 
X [l-i^(X^+2,^)], (A5) 



and 



x^ 



X^M X 



= P(X'^'+",^) X R-'ix^ -r 



X [l-i^(X^4-)] , 
leading to the acceptance factor 



(A6) 



A' 



min < 1 



l-i^:(x^^) i4:(x'^+\ 

P(X'=+2,< 



X(X'=+i,^) 

^('^2M-l) 



(xj+l) l-ii:(Xfe+2,^) 



(A7) 



The generalization to generic s intermediate steps is 
straightforward and can be written as 



A" = min i 1, 



1 - ^^(x^. 



w[x1) 



P(X'=+^^) 

p(xfc,^) 1 - /s:(xfc+^^ 
y:.^ j^(X^+',H ^(^2l/-i) 



(A8) 



To find a simple solution for the unknown probabilities, 
we first impose a cancellation for the intermediate accep- 
tance factors, namely 



i^(X,^) _ wjxi) 



(A9) 



this condition is satisfied by Eqs. and (|20p . Then, 
we notice that the acceptance factor can be written only 
in terms of the final and the initial states as 



A' 



min < 1, 



g(X''~+^^) 



(AlO) 



Further, we can impose the two factors q to be inde- 
pendent on the direction, i.e., the condition (/(X,^) ~ 
q(X, = (/(X), which is satisfied if 



l-if(X,->) 



X wix2M-l) 
P(X,^) 



1 -X(x,^) 



xw(a:i). (All) 



Since the two time directions are mutually exclusive, it is 
also true that P(X,^) -I- P(X,— j') = 1, which allows us 
to solve Eq. (ETT|) and obtain Eqs. (HH) and The 
same reasoning can be repeated for the left direction and, 
due to imposed homogeneity for the probabilities, it can 
be checked that the detailed balance is satisfied for the 
left direction too. 



Appendix B: Bounce algorithm, directed updates, 
and efficiency 

In this Appendix we comment on the relationship be- 
tween the directed-update scheme and the bounce algo- 
rithm. If a = 1 is taken in Eqs. (fT9|) and (|20| . then after 
s updates along the direction d, at the end of the Markov 
step P(X'^+'',(i) — 0, i.e., the next Markov step will be 
taken in the opposite direction, just like the bounce al- 
gorithm. Although the two algorithms are similar in 
this particular limit, there is an important difference 
which eventually leads to a different computational ef- 
ficiency. In order to elucidate this point and to show the 
a-dependence of the efficiency of the directed updates, 
we have done a systematic comparison of the two algo- 
rithms. 

In particular, we have compared the efficiency of the 
directed updates with the bounce algorithm for a one- 
dimensional Heisenberg model. The computational effi- 
ciency is generally defined as 



1 



1 



(Bl) 



where CTq is the square of the statistical error associated 
to a given observable after a given computational time T. 
In Fig. [71 we show the ratio between the directed-update 
scheme efficiency over the bounce algorithm efficiency, 
for the measurement of the ground-state energy of a one- 
dimensional chain. 

We notice that the two sampling schemes have com- 
parable performances, being both based on a similar ap- 
proach. As anticipated, it clearly emerges from Fig. [7] 
that the two algorithms do not have exactly the same 
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1.12 



Figure 7: Relative efficiency of the directed update scfieme 
and the bounce algorithm. The measured quantity is the 
ground-state energy of the one-dimensional Heisenberg model 
on a chain of size L = 80 sites. 



behavior at a = 1, the maximum efficiency of the di- 
rected updates being reached for lower values of a. This 
feature is due to the fact that when a is very close to its 
saturation value, then a single Markov step can consist 
of a conspicuous number of individual "sliding moves". 
Even if this situation leads to a fast decorrelation of con- 
figurations it also leads to a rarefaction of the possibility 
to measure the desired observables, which can eventu- 
ally take place only at the end of the Markov step and 
not during the individual moves. This leads to a worse 
efficiency if compared to the bounce algorithm, where 
measurements can be in principle done after every slid- 
ing move. 

In conclusion, the performances of the two algorithms 
are very close, although some advantages may arise from 
the use of the directed-updates. We further notice that 
the purely Markovian approach introduced in this paper 
could be slightly more efficient in cases where the num- 
ber of rejected configurations by the bounce algorithm is 
substantial whereas all the generated configurations are 
accepted in the directed update scheme. 
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